#Data analysis file for survey experiment carried out in August 2012, examining the effect of information regarding income inequality on attitudes toward income inequality. 
#Kris-Stella Trump kstrump@gmail.com

library(foreign)
#import clean data. During data cleaning process, I removed IP addresses of respondents, removed commas from write-in responses so that these would get read correctly, and replaced "zero" answers on "ought to earn" question (usually given for politicians) with the mean of the remaining responses. 
data <- read.csv("Exp2 - cleandata.csv", na.strings="")
names(data)
dim(data)

#######Creating variables for analysis

#create a variable for republican and democrat, including leaners
#to avoid problems with na, make leaners zero in party.ind only
data$party.ind[is.na(data$party.ind)] <- 0
#indicators of republican, democrat
data$republican <- c()
data$republican <- ifelse((data$party>3 | data$party.ind==3), 1, 0)
data$democrat <- c()
data$democrat <- ifelse((data$party<3 | data$party.ind==1), 1, 0)
#here, both indicators include complete independents in the omitted group

#create treatment variable
data$treat <- (is.na(data$tr.table))==FALSE
sum(data$treat)

#belief in just world
#half the respondents answered this scale before the experimental part of the survey, and half after. this variable combines the pre- and post-bjw scales so that bjw is identified for all participants.
summary(data$bjw.manual) 

#Perception ratios
data$percmin <- apply(data[,c("est_1_1_TEXT", "est_2_1_TEXT", "est_3_1_TEXT", "est_4_1_TEXT", "est_5_1_TEXT", "est_6_1_TEXT")], 1, min)
data$percmax <- apply(data[,c("est_1_1_TEXT", "est_2_1_TEXT", "est_3_1_TEXT", "est_4_1_TEXT", "est_5_1_TEXT", "est_6_1_TEXT")], 1, max)
data$percratio <- (data$percmax/data$percmin)
data$percindex <- log(data$percratio)

#Preference ratios
#allpref are the same type of combination variable as above for bjw. the "allpref" variables include treated and control group respondents. (it is the combination of the pref.noinf and pref.table variables)
data$prefmin <- apply(data[,c("allpref_1", "allpref_2", "allpref_3", "allpref_4", "allpref_5", "allpref_6")], 1, min)
data$prefmax <- apply(data[,c("allpref_1", "allpref_2", "allpref_3", "allpref_4", "allpref_5", "allpref_6")], 1, max)
data$prefratio <- (data$prefmax/data$prefmin)
data$prefindex <- log(data$prefratio)

data$treat <- as.numeric(data$treat)


#### Treatment check: was randomization successful?

#controlling for treatment should not influence perceptions if randomization was successful
t.test(data[data$treat==0,"percratio"], data[data$treat==1,"percratio"])
#no difference between treatment groups in perceptions of inequality


####### Descriptive stats 

dim(data) #407
sum(data$treat) #203

#age (born + 18 for age)
mean(data$born, na.rm=T) #30
median(data$born, na.rm=T) #27
range(data$born, na.rm=T) #19-66

#gender (female=1, male=2)
sum(data$gender==2, na.rm=TRUE) #255 men
sum(data$gender==1, na.rm=TRUE) #149 women

#education - college degree or more
data$degree <- data$educ>3
sum(data$degree, na.rm=T) #223 college degrees

#race 
sum(data$race==1) #314 white
sum(data$race==2) #16 black
sum(data$race>3) #51 other
data$white <- data$race==1

#got through screen on first try
sum(data$att1==5, na.rm=T) #293 

#republican
sum(data$republican, na.rm=T) #80

#bjw
mean(data$bjw.manual) #3.24
range(data$bjw.manual) #1-5.75


##### Balance table

balance <- zelig(treat ~ born + gender + white + republican +degree + bjw.manual + percindex, model="ls", data=data)
summary(balance)

#Age (born+18 for age)
mean(data[data$treat==0,"born"], na.rm=T)
mean(data[data$treat==1,"born"], na.rm=T) 
#Gender
sum(data[data$treat==0 & data$gender==1,"gender"], na.rm=T)/sum(data$treat==0)
sum(data[data$treat==1 & data$gender==1,"gender"], na.rm=T)/sum(data$treat==1)
#Race
sum(data[data$treat==0 & data$white==1,"white"], na.rm=T)/sum(data$treat==0)
sum(data[data$treat==1 & data$white==1,"white"], na.rm=T)/sum(data$treat==1)
#BJW
mean(data[data$treat==0,"bjw.manual"], na.rm=T)
mean(data[data$treat==1,"bjw.manual"], na.rm=T) 
#Republican
sum(data[data$treat==0 & data$republican==1,"republican"], na.rm=T)/sum(data$treat==0)
sum(data[data$treat==1 & data$republican==1,"republican"], na.rm=T)/sum(data$treat==1)
#College degree
sum(data[data$treat==0 & data$degree==1,"degree"], na.rm=T)/sum(data$treat==0)
sum(data[data$treat==1 & data$degree==1,"degree"], na.rm=T)/sum(data$treat==1)
#Perception of inequality
mean(data[data$treat==0,"percindex"], na.rm=T)
mean(data[data$treat==1,"percindex"], na.rm=T)



###### Main results 

#install.packages("Zelig")
library(Zelig)

z.out <- zelig(prefindex ~ treat, model="ls", data=data)
summary(z.out)
#### full regression for comparison of demographic influences 
z.out2 <- zelig(prefindex ~ bjw.manual + treat + republican + percindex, model="ls", data=data)
summary(z.out2)


#Create expected values. 
x.out <- setx(z.out, treat=c(0,1))
s.out <- sim(z.out, x=x.out)
summary(s.out)
#use this output for graphing results


#First graph: raw results, taken from s.out
#mean: treat=0, treat=1, other items at mean
meanb <- c(2.20, 2.67) 
lcib <- c(2.04, 2.50)
ucib <- c(2.36, 2.83)
rawresultsblog <- as.data.frame(cbind(meanb, ucib,lcib))
rawresultsb <- exp(rawresultsblog)

#for error bar command, load hmisc
library(Hmisc)

#First save layout of the plot, to extract midpoints of the bars, used in addition of ci
mp <- barplot(rawresultsb$meanb, space=0.5, ylim=c(0,20), main="Suggested income ratios (highest to lowest salary)", names.arg=c("No information condition", "Information condition"))

pdf("survey1results.pdf")
barplot(rawresultsb$meanb, space=0.5, ylim=c(0,20), main="Recommended income ratios", names.arg=c("No information", "Information condition"), col=c("grey40", "lightgrey"))
axis(1,at=c(0,4),labels=F,lwd.ticks=0)
title(ylab="Recommended income ratio")
errbar(mp,rawresultsb$meanb,rawresultsb$uci,rawresultsb$lci,add=T,pch=NA, )
dev.off()


####Supplemental analysis: comparison of impact by partisan identity
z.out3 <- zelig(prefindex ~ treat*republican, model="ls", data=data)
summary(z.out3)

z.out4 <- zelig(prefindex ~ treat*republican + percindex + bjw.manual, model="ls", data=data)
summary(z.out4)


#### are income differences too large, by treatment

z.out2 <- zelig(opinion_1 ~ treat, model="ls", data=data)
summary(z.out2)

z.out3 <- zelig(opinion_1 ~ treat + republican + bjw.manual + percindex, model="ls", data=data)
summary(z.out3)


####### Additional analysis


#What were average perceptions? 
summary(data$percratio)
#actual gap is
11400000/24240 #470.297
sum(data$percratio>470) #29
29/407 #0.07, ie 7% are right or above in estimate

#Average prefs without treatment
control <- data[data$treat==0,] #204
sum(control$prefratio>=470) #1
1/407 #0.2% ok with current division
summary(control$prefratio)
#3rd quartile is 30
quantile(control$prefratio, probs=c(0.025,0.5,0.95,0.975))
#95% want a ratio of 65 or below

#Prefs after treatment
treatment <- data[data$treat==1,]
sum(treatment$prefratio>=470) #2 people ok with current division
summary(treatment$prefratio)


#Treatment condition and opinion questions

#does treatment influence trust
t.test(data[c(data$treat==1), "trust.gov"], data[c(data$treat==0), "trust.gov"])
#p=0.44

#does treatmet influence assessment of income diff too large
t.test(data[c(data$treat==1), "opinion_1"], data[c(data$treat==0), "opinion_1"])
#p=0.13, effect size 0.15

#large diff necessary
t.test(data[c(data$treat==1), "opinion_2"], data[c(data$treat==0), "opinion_2"])
#p=0.20, effect size 0.13

#resp of gov to reduce income gap
t.test(data[c(data$treat==1), "opinion_3"], data[c(data$treat==0), "opinion_3"])
#p=0.14, effect size 0.18

#rich pay too much in taxes
t.test(data[c(data$treat==1), "opinion_4"], data[c(data$treat==0), "opinion_4"])
#p=0.55

#gov resp to help the poor
t.test(data[c(data$treat==1), "opinion_5"], data[c(data$treat==0), "opinion_5"])
#p=0.54

